home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / vol_200 / 275_02 / prob21.c < prev    next >
Text File  |  1980-01-01  |  29KB  |  1,018 lines
  1.  
  2. /* prob21.c                            */
  3. /* program for lcau21 option 't'                */
  4. /* calculate probabilities related to evolution            */
  5. /* Harold V. McIntosh, 10 August 1987                */
  6.  
  7. /* references:                            */
  8. /*                                */
  9. /* W. John Wilbur, David J. Lipman and Shihab A. Shamma        */
  10. /* On the prediction of local patterns in cellular automata    */
  11. /* Physica 19D 397-410 (1986)                    */ 
  12. /*                                */
  13. /* Howard A. Gutowitz, Jonathan D. Victor and Bruce W. Knight    */
  14. /* Local structure theory for cellular automata            */
  15. /* Physica 28D 18-48 (1987)                    */
  16.  
  17. /*    Copyright (C) 1987    */
  18. /*    Copyright (C) 1988    */
  19. /*    Harold V. McIntosh    */
  20. /*    Gerardo Cisneros S.    */
  21.  
  22. # define BROW        3    /* row for bar charts        */
  23. # define EROW        1    /* row for evolution synopsis    */
  24. # define AORG         0    /* x-origin  of contour plot    */
  25. # define BORG       109    /* x-origin  of 2-block param   */
  26. # define CORG       219    /* x-origin  of Bernstein plot  */
  27.  
  28. /* edit the probability screen    */
  29. edtri() {char c;
  30.   videomode(COLGRAF);
  31.   videopalette(YELREGR);
  32.  
  33.   while (0<1) {
  34.   woruno(0,28);
  35.   videocursor(0,0,36);
  36.   videoputc('?',2);
  37.   c=kbdin();
  38.   if (c == '\015') break;
  39.   videocursor(0,0,38);
  40.   videoputc(c,2);
  41.   videocursor(0,0,36);
  42.   videoputc(' ',0);
  43.   switch (c) {
  44.     case '+': videopalette(WHCYMAG); break;
  45.     case '-': videopalette(YELREGR); break;
  46.     case 'a': asfreq(3); break;
  47.     case 'e': pevolve(); break;
  48.     case 'g': lifreq(50,2); break;
  49.     case 'G': lifreq(200,1); break;
  50.     case 'm': moncar(1,2); break;
  51.     case 'M': moncar(50,1); break;
  52.     case 'x': pdiff(100); break;
  53.     case 'i': pidiff(100); break;
  54.     case 'j': pjdiff(100); break;
  55.     case 'y': pdiff3(100); break;
  56.     case 'z': pdiff4(100); break;
  57.     case 'w': pdiff5(100); break;
  58.     case 'v': pdiff6(100); break;
  59.     case 't': twoblockfreq(100); break;
  60.     case '1': nblclr(); oneblfreq(8*BROW,300,48); break;
  61.     case '2': nblclr(); twoblfreq(8*BROW,300,48); break;
  62.     case '3': nblclr(); thrblfreq(8*BROW,300,48); break;
  63.     case '4': nblclr(); foublfreq(8*BROW,300,48); break;
  64.     case '5': nblclr(); fivblfreq(8*BROW,300,48); break;
  65.     case '6': nblclr(); sixblfreq(8*BROW,300,48); break;
  66.     case '/': videomode(COLGRAF); videopalette(YELREGR); break;
  67.     case '?': trmenu(); break;
  68.     default: break;
  69.     }; /* end switch */
  70.   };   /* end while  */
  71.   videopalette(WHCYMAG);
  72.   videomode(T80X25);
  73. }      /* end edtri  */
  74.  
  75. /* show menu */
  76. trmenu() {
  77.   videoscroll(BROW,0,BROW+8,40,0,0);
  78.   videocursor(0,BROW,0);
  79.   printf("a       - a priori estimates\n");
  80.   printf("m,M,g,G - sample evolution\n");
  81.   printf("xyzwv   - selfconsistent probabilities\n");
  82.   printf("xij     - iterated s-c probabilities\n");
  83.   printf("t       - graph 2 block probabilities\n");
  84.   printf("123456  - n-block bar charts\n");
  85.   printf("+-      - change color pallette\n");
  86.   printf("e       - 12 lines evolution\n");
  87.   printf("/?(clear screen, show menu), <cr>(exit)\n");
  88. }
  89.  
  90. /* show fourteen lines of evolution at top of screen */
  91. pevolve() {int i, j;
  92.   videoscroll(EROW,0,EROW+1,40,0,0);
  93.   asctobin();
  94.   for (j=8*EROW; j<8*(EROW+2)-2; j++) {
  95.     for (i=0; i<AL; i++) videodot(j,i,arr1[i]);
  96.     onegen(AL);
  97.     };
  98. }
  99.  
  100. /* Clear a space for the n-block frequencies */
  101. nblclr() {videoscroll(BROW,0,BROW+8,40,0,0);}
  102.  
  103. /* make a frame for a graph          */
  104. /* (x,y) = lower left corner; e.g. (0,0) */
  105. /* n     = dimension of frame            */
  106. gfram(x,y,n) int x, y, n; {int i;
  107.  
  108. for (i=0; i<=n; i++) videodot(199-y-i,x,0);
  109. for (i=0; i<=n; i++) videodot(199-y-i,x+n,0);
  110. for (i=0; i<=n; i++) videodot(199-n-y,x+i,0);
  111. for (i=0; i<=n; i++) videodot(199-y,x+i,0);
  112.  
  113. for (i=0; i<=n; i+=10) videodot(199-y-i,x,3);
  114. for (i=0; i<=n; i+=10) videodot(199-y-i,x+n,3);
  115. for (i=0; i<=n; i+=10) videodot(199-n-y,x+i,3);
  116. for (i=0; i<=n; i+=10) videodot(199-y,x+i,3);
  117. }
  118.  
  119. /* put a diagonal in a graph */
  120. gdiag(x,y,n) int x, y, n; {int i;
  121. for (i=0; 2*i<n; i++) videodot(199-y-2*i,x+2*i,3);
  122. }
  123.  
  124. /* graph Bernstein polynomial */
  125. bgraf(x,y,k,n) int x, y, k, n; {int i; double bern(), en, p;
  126. en=(double)(n);
  127. for (i=0; i<n; i++) {
  128.   p=((double)(i))/en;
  129.   videodot(199-y-(int)(en*bern(p,k)),x+i,1);
  130.   };
  131. }
  132.  
  133. /* "Monte Carlo" determination of probabilities */
  134. moncar(n,l) int n, l; {
  135. int i, j, k, b[KK], bb[KK][KK];
  136. double bf[KK], bbf[KK][KK];
  137.  
  138. nblclr();
  139. gfram(BORG,0,100);
  140. asctobin();
  141. for (k=0; k<n; k++) {
  142.   onegen(AL);
  143.   for (i=0; i<KK; i++) b[i]=0;
  144.   for (i=0; i<AL; i++) b[arr1[i]]++;
  145.   for (i=0; i<KK; i++) bf[i]=((double)(b[i]))/((double)(AL));
  146.   for (i=0; i<KK; i++) for (j=0; j<KK; j++) bb[i][j]=0;
  147.   for (i=1; i<AL; i++) bb[arr1[i-1]][arr1[i]]++;
  148.   bb[arr1[AL-1]][arr1[0]]++; 
  149.   for (i=0; i<KK; i++) for (j=0; j<KK; j++) 
  150.     bbf[i][j]=((double)(bb[i][j]))/((double)(AL));
  151.   videodot(199-(int)(100.0*bbf[1][1]),BORG+(int)(100.0*bf[1]),l);
  152.   };
  153.  videocursor(0,BROW+7,0);
  154. printf("(Monte Carlo) "); 
  155. for (i=0; i<KK; i++) printf("%2d:%5.3f ",i,bf[i]);
  156. videocursor(0,BROW+8,0);
  157. for (i=0; i<KK; i++) for (j=0; j<KK; j++) 
  158.   printf("%1d%1d:%5.3f ",i,j,bbf[i][j]);
  159. }
  160.  
  161. /* Generate coefficients of 2nd generation Bernstein Polynomial */
  162. berncoef() {
  163. int i, i0, i1, i2;
  164.  
  165. for (i=0; i<BD; i++) bp[i]=0.0;
  166. for (i0=0; i0<KK; i0++) {
  167. for (i1=0; i1<KK; i1++) {
  168. for (i2=0; i2<KK; i2++) {
  169.   if (ascrule[i0][i1][i2]=='1') bp[i0+i1+i2]+=1.0;
  170.   };};};
  171. }
  172.  
  173. /* Generate coefficients of 3rd generation Bernstein Polynomial */
  174. bernthrd() {
  175. int i, i0, i1, i2, i3, i4;
  176. int j0, j1, j2;
  177.  
  178. for (i=0; i<BD; i++) bp[i]=0.0;
  179. for (i0=0; i0<KK; i0++) {
  180. for (i1=0; i1<KK; i1++) {
  181. for (i2=0; i2<KK; i2++) {
  182. for (i3=0; i3<KK; i3++) {
  183. for (i4=0; i4<KK; i4++) {
  184.   j0=ascrule[i0][i1][i2]-'0';
  185.   j1=ascrule[i1][i2][i3]-'0';
  186.   j2=ascrule[i2][i3][i4]-'0';
  187.   if (ascrule[j0][j1][j2]=='1') bp[i0+i1+i2+i3+i4]+=1.0;
  188.   };};};};};
  189. }
  190.  
  191. /* Generate coefficients of 4th generation Bernstein Polynomial */
  192. bernfrth() {
  193. int i, i0, i1, i2, i3, i4, i5, i6;
  194. int j0, j1, j2, j3, j4;
  195. int k0, k1, k2;
  196.  
  197. for (i=0; i<BD; i++) bp[i]=0.0;
  198. for (i0=0; i0<KK; i0++) {
  199. for (i1=0; i1<KK; i1++) {
  200. for (i2=0; i2<KK; i2++) {
  201. for (i3=0; i3<KK; i3++) {
  202. for (i4=0; i4<KK; i4++) {
  203. for (i5=0; i5<KK; i5++) {
  204. for (i6=0; i6<KK; i6++) {
  205.   j0=ascrule[i0][i1][i2]-'0';
  206.   j1=ascrule[i1][i2][i3]-'0';
  207.   j2=ascrule[i2][i3][i4]-'0';
  208.   j3=ascrule[i3][i4][i5]-'0';
  209.   j4=ascrule[i4][i5][i6]-'0';
  210.   k0=ascrule[j0][j1][j2]-'0';
  211.   k1=ascrule[j1][j2][j3]-'0';
  212.   k2=ascrule[j2][j3][j4]-'0';
  213.   if (ascrule[k0][k1][k2]=='1') bp[i0+i1+i2+i3+i4+i5+i6]+=1.0;
  214.   };};};};};};};
  215. }
  216.  
  217. /* Generate coefficients of 5th generation Bernstein Polynomial */
  218. bernfifth() {
  219. int i, i0, i1, i2, i3, i4, i5, i6, i7, i8;
  220. int j0, j1, j2, j3, j4, j5, j6;
  221. int k0, k1, k2, k3, k4;
  222. int l0, l1, l2;
  223.  
  224. for (i=0; i<BD; i++) bp[i]=0.0;
  225. for (i0=0; i0<KK; i0++) {
  226. for (i1=0; i1<KK; i1++) {
  227. for (i2=0; i2<KK; i2++) {
  228. for (i3=0; i3<KK; i3++) {
  229. for (i4=0; i4<KK; i4++) {
  230. for (i5=0; i5<KK; i5++) {
  231. for (i6=0; i6<KK; i6++) {
  232. for (i7=0; i7<KK; i7++) {
  233. for (i8=0; i8<KK; i8++) {
  234.   j0=ascrule[i0][i1][i2]-'0';
  235.   j1=ascrule[i1][i2][i3]-'0';
  236.   j2=ascrule[i2][i3][i4]-'0';
  237.   j3=ascrule[i3][i4][i5]-'0';
  238.   j4=ascrule[i4][i5][i6]-'0';
  239.   j5=ascrule[i5][i6][i7]-'0';
  240.   j6=ascrule[i6][i7][i8]-'0';
  241.   k0=ascrule[j0][j1][j2]-'0';
  242.   k1=ascrule[j1][j2][j3]-'0';
  243.   k2=ascrule[j2][j3][j4]-'0';
  244.   k3=ascrule[j3][j4][j5]-'0';
  245.   k4=ascrule[j4][j5][j6]-'0';
  246.   l0=ascrule[k0][k1][k2]-'0';
  247.   l1=ascrule[k1][k2][k3]-'0';
  248.   l2=ascrule[k2][k3][k4]-'0';
  249.   if (ascrule[l0][l1][l2]=='1') bp[i0+i1+i2+i3+i4+i5+i6+i7+i8]+=1.0;
  250.   };};};};};};};};};
  251. }
  252.  
  253. /* Generate coefficients of 6th generation Bernstein Polynomial */
  254. bernsixth() {
  255. int i, i0, i1, i2, i3, i4, i5, i6, i7, i8, i9, i10;
  256. int j0, j1, j2, j3, j4, j5, j6, j7, j8;
  257. int k0, k1, k2, k3, k4, k5, k6;
  258. int l0, l1, l2, l3, l4;
  259. int m0, m1, m2;
  260.  
  261. for (i=0; i<BD; i++) bp[i]=0.0;
  262. for (i0=0; i0<KK; i0++) {
  263. for (i1=0; i1<KK; i1++) {
  264. for (i2=0; i2<KK; i2++) {
  265. for (i3=0; i3<KK; i3++) {
  266. for (i4=0; i4<KK; i4++) {
  267. for (i5=0; i5<KK; i5++) {
  268. for (i6=0; i6<KK; i6++) {
  269. for (i7=0; i7<KK; i7++) {
  270. for (i8=0; i8<KK; i8++) {
  271. for (i9=0; i9<KK; i9++) {
  272. for (i10=0; i10<KK; i10++) {
  273.